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1. Introduction 

Stars are hot plasmas. Therefore, the presence of a solid neutron star (NS) crust may seem 
unlikely. Nevertheless compact stars, such as white dwarfs (WD) and NS, are so dense that 
this plasma can crystallize. How do these stars freeze, and what are some of the extraor- 
dinary properties of this solid star stuff? In this chapter we describe molecular dynamics 
(MD) simulations of plasma crystals and NS crust. These simulations determine many crust 
properties that may be crucial for the electromagnetic, neutrino, and gravitational wave ra- 
diations of NS. 
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We start in Sec. [T] describing how plasmas freeze in the laboratory, in the interior of 
white dwarf stars, and to form new neutron star crust. In Sec. |3]] we describe a molecular 
dynamics formalism for accurately determining many crust properties. In Sec. |4|] we review 
MD results for chemical separation upon freezing, diffusion in Coulomb solids, breaking 
strain (strength of the NS crust), shear viscosity, and dynamical response. Finally, Sec. [5] 
contains a short summary and discusses some open questions and future challenges. 

2. How do Stars Freeze? 

How do compact stars freeze? Surprisingly there has been little work on this fundamental 
question. In this section we first examine crystallization in laboratory plasmas and then 
discuss freezing in white dwarfs and neutron stars. 

2.1. Crystallization of Laboratory Plasmas 

In the laboratory, one can study complex (or dusty) plasma crystals. Complex plasmas (CP) 
are low temperature plasmas containing charged microparticles, for a review see Fortov et 
al. fT] . Often the microparticles are micron sized spheres that acquire large electric charges 
and the strong coulomb interactions between microparticles can lead to crystallization. In- 
deed plasma crystals were first observed in the laboratory in 1994 [2]. Complex plasmas 
typically differ from White Dwarf interiors and Neutron Star crusts in a number of ways. 
First the microparticles feel additional fluctuating and friction forces because of interactions 
with the background gas. Note that in stars, electron-ion interactions are small because of 
the large electron degeneracy. Second, the Debye screening length A, see Eq. [T]in Sec. [3jJ 
is often smaller in the CP than in a star (when measured in units of the lattice spacing). 
This changes the lattice type from body-centered-cubic (bcc) as expected in stars, to face- 
centered-cubic (fee) or other types in a CP. Finally in a CP there is an overall confining 
potential, and because of gravitational gradients it is often easier to study two-dimensional 
CP crystals. 

Three-dimensional CP crystals have been formed onboard the International Space Sta- 
tion under microgravity conditions. Details of the experiment are presented in ref. [3|. 
Microparticles were found in regions with fee and hexagonal-close packing (hep) order 
in 21, in agreement with MD simulations [4]. Khrapak et al. |[6l studied freezing and melt- 
ing of these CP crystals and found diffusion to be relatively fast so that the system remained 
in equilibrium. Melting criteria for CP systems were presented by Klumov [7J . It should 
be possible to record complete trajectories, position as a function of time, for each of the 
microparticles in a three dimensional CP experiment. This would allow a very detailed 
dynamical study of the melting or freezing phase transitions. 

2.2. Crystallization in White Dwarf stars 

Observations of cooling White Dwarf (WD) stars provide important information on the ages 
and evolution of stellar systems ||8l|9l|T0lllT|. The interior of a WD is a coulomb plasma 
of ions and a degenerate electron gas. As the star cools this plasma crystallizes. Note 
that WDs freeze from the center outward, because of the higher central densities. While 
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NS are expected to form a solid crust over a liquid core. This crystallization can delay 
WD cooling as the latent heat is radiated away, see for example ref. |[T2ll . Winget et al. 
recently observed effects from the latent heat of crystallization on the luminosity function 
of WDs in the globular cluster Ngc 6397 [13]. Winget et al.'s observations may constrain 
the melting temperature of the carbon and oxygen mixtures expected in these WD cores. In 
addition, Astro-seismology provides an alternative way to study crystallization in WD, see 
for example lfT4l . 

How do stars freeze? At present we have very little direct observational information to 
answer this fundamental question. In the future we can observe many more WD systems 
and this should greatly increase the observational information on crystallization in WD. 



2.3. Crystallization on Accreting Neutron Stars 

We next consider freezing in accreting NS. Material falling on a NS can undergo rapid pro- 
ton capture (or rp process) nucleosynthesis to produce a range of nuclei with mass numbers 
A that could be as high as A 100 lITSl [161 . As this rp process ash is buried by further 
accretion, the rising electron fermi energy induces electron capture to produce a range of 
neutron rich nuclei from O to approximately Se lITTll . This material freezes, when the den- 



sity reaches near 10^*^ g/cm^. In Sec. 4.1. we describe large scale MD simulations of how 
this complex rp process ash freezes [18|. We find that chemical separation takes place and 
the liquid ocean is greatly enriched in low atomic number Z elements, while the newly 
formed solid crust is enriched in high Z elements. Note that for an accreting star, there is 
typically a thin liquid ocean atop a solid crust that itself surrounds a dense liquid core. 

Furthermore, MD simulations, as discussed in Sec. |4]] suggest that a regular crystal 
lattice forms even though large numbers of impurities may be present. This regular crystal 
should have a high thermal conductivity. We do not find an amorphous solid that would 
have a low thermal conductivity lfT9l . 

Recent X-ray observations of neutron stars, find that the crust cools quickly, when heat- 
ing from extended periods of accretion stops |[20l 1211 l22l . This is consistent with MD 
simulations, and strongly favors a crystalline crust over an amorphous solid that would cool 
more slowly lf23l. Il24]|. 



3. Molecular Dynamics Formalism 

Neutron star crust consists of a relativistic fermi gas of electrons, a crystal lattice of neutron 
rich ions, and, in general, a neutron gas. For a review see |[25l . In Section [3X] we de- 
scribe our MD simulation formalism and then discuss some computational considerations 
in Section [3^ 



3.1. Classical Molecular Dynamics for Coulomb plasmas 

We model NS crust as a classical Coulomb plasma. The potential between the i*^ and j*^ 
ions is assumed to be of screened Coulomb (Yukawa) form, 

v^,{r) = ^e-A. (1) 



4 



C. J. Horowitz et al. 



Here the ions are separated by a distance r and have charges Zi and Zj. The Fermi screening 
length A, for cold relativistic electrons, is A^^ = 2a^/'^kp /tt^/"^ where the electron Fermi 
momentum kp is kp = (Svr^ne)^/^ and a is the fine structure constant. The electron 
density rig is equal to the ion charge density, ng = {Z)n, where n is the ion density and 
{Z) is the average charge. Our simulations are classical and we have neglected the electron 
mass (extreme relativistic limit). This is to be consistent with our previous work on neutron 
stars. However, the electron mass is important at the lower densities in WD and this may 
change our results slightly [26]. We assume the ions are classical. In general the large mass 
of the ions ensures that their thermal wavelengths are smaller than the inter-particle spacing. 
However, quantum effects could play some role at high densities EtI . ESII . For relativistic 
electrons, the ratio of A to the ion sphere radius a, 

/ 3 n1/3 

depends only on the average charge {Z). 

The simulations can be characterized by a Coulomb parameter T describing the ratio 
of a typical Coulomb energy to the thermal energy. For a one component plasma (OCP), 
where all of the ions have charge Z, at a temperature T we have, 

Z\^ 

For a mixture of ions of different charges we define F using an appropriate average of Zf / ai 
where the typical distance between ions ai is expected to scale as z]^^ so that, 

r = ^ ^ . (4) 

Here {Z^/^) is an average over the ion charges, T is the temperature, and the electron sphere 
radius Og is ag = (3/47rng)^/'^ with rig = {Z)n the electron density. The one component 
system freezes near T = 175, while mixtures of ions are expected to freeze at a somewhat 
higher r imill. 

Time can be measured in units of one over the plasma frequency ojp. Long wavelength 
fluctuations in the charge density can undergo oscillations at the plasma frequency. This 
depends on the ion charge Z and mass M. For mixtures we define a hydrodynamical 
plasma frequency Ojp from the simple averages of Z and M, 

r47re2(Z)2nii/2 



3.2. Computational Considerations for MD Simulations 

The pure unscreened Coulomb interaction has an effectively infinite range. MD simulations 
involving pure Coulomb interactions must take this into account to produce physically ac- 
curate results. Often this is done by assuming periodic boundary conditions and using 
Ewald summation to account for interactions with distant particles. At the opposite ex- 
treme, systems with very short-range interactions, such as the Lennard- Jones potential, can 
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be simulated with an MD code that uses a very small cut-off distance. One assumes the 
interaction is zero for distances beyond the cutoff. In these systems, each particle interacts 
with only a few tens to hundreds of nearby neighbors, even if the system as a whole in- 
volves millions of particles. The screened Coulomb interaction in our ion simulations falls 
between these extremes. We do not need to account for interactions with very distant ions, 
yet because the screening length A is rather large, a small cut-off radius will not model the 
physics properly. We find that a cutoff of about 8A gives good physical results. Note that a 
cutoff distance significantly smaller than 8A will lead to errors in the shear modulus |[3T| . 
We use a velocity Verlet integration scheme to advance the system in time [32|. 

Our computer code is a parallel Fortran 95 program that uses both MPI (Message Pass- 
ing Interface) and OpenMP. The force calculation is carried out in a pair of nested DO 
loops. The outer loop is over what we call target particles, while the inner loop is over 
source particles. These are the same sets of particles of course, but for sake of discussion it 
is convenient to distinguish the particles that are being acted on and those doing the acting. 

Most supercomputers at the time of this writing have a hybrid architecture. These ma- 
chines consist of many individual compute nodes which are complete computers within 
themselves. Each node in turn consists of two or more processor chips sharing a common 
physical memory. Processors chips in turn have several cores that share a common Level 
2 cache and interface to main memory. Each core has its own program counter, floating 
point and integer registers, and floating point, integer and logical units. Each core is there- 
fore capable of running a thread, or light-weight process. This type of node architecture is 
similar to a symmetric multiprocessor (SMP), the kind of architecture for which OpenMP 
was designed. On the other hand, MPI is designed for multi-node machines, where data can 
be shared among nodes only by passing messages. The combination of MPI and OpenMP 
in one code is often the best programming paradigm for multi-node systems where each 
node is an SMP. For example, we have run our MPI-i-OpenMP code on a Cray XT5 com- 
puter called Kraken, which consists of 9408 or more nodes, each of which has two six-core 
processors. We typically use 64 nodes, placing an MPI process on each of node's two pro- 
cessors, and assigning on OpenMP thread to each processor's six cores. This makes a total 
of 768 threads. Each MPI process gets a complete set of particles. Targets are distributed 
among MPI processes, while sources are distributed among OpenMP threads. Each process 
lets all sources act on its targets. 

Our code scales remarkably well. Runs using 768 cores of Kraken obtain almost linear 
speedup on the 55296-ion runs, and only slightly worse on 27648-ion runs. Our longest 
simulation times, of order 2 x 10^/ujp, involved a week or two of computer time with 768 



cores. Note that the breaking strain simulations in Sec. 4.3. used a modified version of the 
program SPaSM 133J. 

4. MD Simulation Results 

In this section we review MD simulation results for a number of neutron star crust proper- 



ties. We start in Sec. 4.1. with chemical separation that occurs as new crust forms on an 



accreting neutron star. Next, in Sec. 4.2. we describe diffusion in Coulomb crystals that 



may be important for the structure of NS crust. The breaking strain (strength) of NS crust 



is discussed in Sec. 4.3. This may be important for mountains on rotating NS that radiate 
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gravitational waves. The shear viscosity of nuclear pasta is discussed in Sec. |4.4.| This 



may be important for the damping of collective modes. Finally we discuss the dynamical 



response of nuclear pasta in Sec. 4.5 



4.1. Phase separation in the crust of accreting neutron stars 

Phase separation is important for white dwarf stars ||29]| . As a star cools, and crystallization 
takes place, the crystal phase is enriched in oxygen while the liquid is enriched in carbon. 
Phase separation is also important for neutron stars that accrete material from a companion. 
This material can undergo nuclear reactions involving rapid proton capture (the rp process) 
to synthesize a variety of medium mass nuclei |[T5]| . Further accretion increases the den- 
sity of a fluid element until crystallization occurs near a density of 10^*^ g/cm^. However, 
molecular dynamics simulations show that this crystallization is accompanied by chemical 
separation. The composition of the new solid crust is very different from the remaining 
liquid ocean. This changes many properties of the crust and can impact many observables. 




Figure 1. Final configuration of carbon ions (light red) and oxygen ions (black) in a 55296 
ion simulation that consisted of 75% oxygen and 25% carbon |[30l . Image prepared with 
VMD 0411 . 

We start with the simpler two component carbon-oxygen system in the interior of a 
white dwarf ll29l [30l . To determine the melting temperature and the composition of the 
liquid and solid phases that are in equilibrium, we perform two-phase MD simulations. A 
region of solid phase is combined with a region of liquid phase, in a rectangular simulation 
volume, as shown in Fig. [T] Here a regular solid phase is on the left and the liquid phase is 
on the right. With periodic boundary conditions, there are two liquid-solid interfaces, one 
near the left side and one near the center. 

This system is then evolved for a long time while the temperature is continually ad- 
justed to keep about half of the system solid and half liquid. During this time carbon ions 
can diffuse into or out of the soUd until the two phases reach chemical equilibrium. The 
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composition of the liquid and solid phases that are in equilibrium can be determined by sim- 
ply counting the number of ions that remain in different regions of the simulation volume. 
The resulting phase diagram is shown in Fig. [2] The melting temperature of an approxi- 
mately 50% carbon, 50% oxygen mixture, that is expected in WD interiors, is seen to be 
only slightly higher than the melting temperature of pure carbon. This result appears to be 
in good agreement with observations of WD cooling in the globular cluster NGC6397 [,13 j . 
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Figure 2. Carbon-oxygen phase diagram plotting the composition of the liquid phase (upper 
red curve or circles) that is in equilibrium with the solid phase (lower blue curve or squares) 
at a melting temperature T in units of the melting temperature for pure carbon Tq, from 
ref. QUI . The number fraction of oxygen ions is xq- Results from 55296 ion simulations 
are filled symbols while the open symbols are results with 27648 ions. The curves are from 
the analytic model of Medin and Cumming ll35l . 

We now consider liquid-solid phase equilibria and chemical separation for more com- 
plicated compositions that may be present as accreting neutron stars form new crust. With 
chemical separation, the liquid ocean is greatly enriched in low atomic number Z elements, 
see Fig. [3] Carbon, if present, may be depleted in the crystal (crust) and enriched in the 
liquid ocean phase. Also, chemical separation may change the thermal conductivity of 
the crust and its temperature profile. Indeed some neutron stars are observed to produce 
energetic X-ray bursts known as superbursts. These are thought to involve the unstable 
thermonuclear burning of carbon Il36ll37l[38l . However it is unclear how the initial carbon 
concentration is obtained and how the ignition temperature is reached. 

Schatz et al. have calculated the rapid proton capture (rp) process of hydrogen burning 
on the surface of an accreting neutron star [15j . This produces a variety of nuclei up to 
atomic masses A ICQ. Gupta et al. |[39l then calculate how the composition of the rp 
process ash evolves because of electron capture and light particle reactions as the material 
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Figure 3. Schematic diagram of the surface of an accreting neutron star fW\. Chemical 
separation upon crystaUization will take place at the boundary between the liquid ocean 
and the solid crust. The ocean is enriched in low Z elements. 



is buried by further accretion. Their final composition, at a density of 2.16 x 10^^ glcw? 
(near neutron drip at the bottom of the outer crust) has forty percent of the ions with atomic 
number Z = 34, while an additional 10 % have Z = 33. The remaining half of the ions 
have a range of lower Z from Z = 8 to 32. Finally there is a small abundance of Z = 36 
and Z = 47. This complex composition is shown in Fig. |4] 

We have performed a two phase MD simulation using this complex 17 component rp- 
ash composition ifTSl . The MD simulation is similar to the above carbon-oxygen simulation 
except that only 27648 ions were used. The results are shown in Fig. [4] The liquid ocean is 
seen to be greatly enriched in low Z elements such as oxygen compared to the solid crust. 
Medin and Gumming have shown how this chemical seperation upon crystallization can 
lead to enrichment of the entire ocean HOll . 

4.2. Diffusion 

Diffusion of impurities, dislocations, or other defects can be very important for the structure 
of neutron star crust and could impact transport and mechanical properties. Diffusion in 
coulomb plasma liquids has been well studied 11411 and is important for sedimentation of 
impurities in white dwarf (WD) ||42l|43llBl and neutron stars (NS) f45^. '46 |. Here ions, 
with a larger than average mass to charge ratio, sink in a strong gravitational field. This 
releases gravitational energy that can delay the cooling of metal rich WD BTl . However, 
we are not aware of previous numerical results for diffusion constants of coulomb crystals 
under Astrophysical conditions. Often the diffusion constant is simply assumed to be zero. 
This diffusion could be important for sedimentation in soUd WD interiors, over long time 
scales, and for the structure of NS crusts. 

Solid diffusion can depend dramatically on the form of the interaction between parti- 
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Figure 4. Abundance (by number) of chemical elements versus atomic number Z for rp 
process ash material on an accreting NS, from ref. lilSJ . The plus symbols show the initial 
composition of the mixture. The final compositions of the liquid phase, open green circles, 
and solid phase, filled red squares, are also shown. 



cles and may be very slow for hard-core systems. For example, the binary Lennard Jones 
(LJ) system with a hard-core oc r~^^ interaction forms a glass because of very slow diffu- 
sion P8l . In contrast, the coulomb plasma with a soft 1/r core should have much faster 
diffusion. Therefore the Coulomb crystal may provide an important model system where 
diffusion is fast enough to be more easily studied by molecular dynamics (MD) simulations. 

We have studied diffusion in solid one component plasmas [19 j. For example. Fig. 
[5] shows ions that have diffused in a small 3456 ion simulation that started from a nearly 
perfect body centered cubic lattice. The diffusing ions are seen to form a chain where 
ions hop to fill vacancies formed by other diffusing ions. Neutron star crust is under great 
pressure and this suppresses vacancy formation. Therefore ions move in ways to minimize 
vacancies. 

Next in Fig. |6] we show diffusion in a 27648 ion system consisting of two micro- 
crystals of different orientations. Ions are seen to diffuse quickly along the grain boundaries 
between micro-crystals. Finally we considered diffusion in amorphous systems where a 
liquid configuration was quickly quenched from the melting temperature to a much lower 
temperature. We find that diffusion in these systems is fast enough so that the systems can 
start to crystalize even over the relatively short time scales accessible to MD simulations. 
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Figure 5. Diffusion in a sample configuration of 3456 ions at F = 175 |[T9l . Ions that have 
moved less than 1.34a in a time t = 236/a;p are small brown dots. Ions that have moved 
more than 1.34a are shown as larger black disks and are seen to be in a ring configuration 
where ions "hop" to lattice sites vacated by other hopping ions. This system started from a 
perfect bcc lattice. Figure plotted with VMD [34.1 . 



This strongly suggests that coulomb solids in WD interiors and NS crust will be crystalline 
and not amorphous. 

In summary, diffusion in Coulomb crystals is relatively fast because the ions have soft 
1/r cores and can slide past one another. As a result astrophysical solids in WD and NS are 
likely to be nearly perfect crystals with very few defects. This suggests they should have 
a high thermal conductivity. This is consistent with observations of rapid crust cooling of 
neutron stars following extended periods of accretion 11211 122ll23ll24l . This rapid cooling 
implies a high crust thermal conductivity, that agrees with the conductivity of a regular 
crystal, and is larger than the conductivity expected for an amorphous solid. 

4.3. Breaking Strain and Gravitational Waves 

The breaking strain (strength) of NS crust determines the maximum sized mountain that is 
possible on a NS before it collapses under the stars extreme gravity. A large mountain, on a 
rapidly rotating NS, produces a time dependent mass quadrupole moment. This efficiently 
radiates gravitational waves. Albert Einstein, almost 100 years ago, predicted the oscillation 
of space and time known as gravitational waves (GW). Within a few years, with the oper- 
ation of Advanced LIGO [49], Advanced VIRGO [.50,1 and other sensitive interferometers, 
we anticipate the historic detection of GW. 

The first GW that are detected will likely come from the merger of two neutron stars. 
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Figure 6. Diffusion in a configuration of 27648 ions starting from imperfect crystal initial 
conditions lfT9l . Ions that move only a small distance are small gray points. Ions that have 
moved over three lattice spacings, during the simulation time of t = SOOOO/wp, are shown 
as large blue spheres. These are seen to be clustered at the grain boundaries. The initial 
conditions included two micro-crystals of different orientation. Figure plotted using VMD 



The rate of such mergers can be estimated from known binary systems lISTll . During a 
merger the GW signal has a so called chirp form where the frequency rises as the two 
neutron stars spiral closer together. Deviations of this wave form from that expected for 
two point masses may allow one to deduce the equation of state of neutron rich matter and 
measure the radius of a neutron star r^s ll52l . Alternatively one may be able to observe 
the frequency of oscillations of the hyper-massive neutron star just before it collapses to 
a black hole. This frequency depends on the radius of the maximum mass neutron star 
|[53l . However, either approach may require high signal to noise data from relatively nearby 
mergers. 

Continuous GW signals can also be detected, see for example ||54]| . in addition to burst 
signals such as those from neutron star mergers. There are several very active ongoing 
and near future searches for continuous gravitational waves at LIGO, VIRGO and other 
detectors, see for example [55 1. Often one searches at twice the frequency of known ra- 
dio signals from pulsars because of the quadrupole nature of GW. No signal has yet been 
detected. However, sensitive upper limits have been set. These limits constrain the shape 
of neutron stars. In some cases the star's elipticity e, which is that fractional difference in 
moments of inertia e = (/i — h)/^ is observed to be less than a part per million or even 
smaller. Here h, h, and I3 are the principle moments of inertia. 

In general, the amplitude of any continuous signal is much weaker than a burst signal. 
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However, one can gain sensitivity to a continuous signal by coherently (or semi-coherently) 
integrating over a large observation time, see for example ll56l . Note that searches for 
continuous GW can be very computationally intensive because one must search over an 
extremely large space of parameters that may include the source frequency, how that fre- 
quency changes with time, the source location on the sky, etc. The Einstein at home dis- 
tributed computing project uses spare cycles on the computers of a large number of volun- 
teers to search for continuous GW lISTll . 

Strong GW sources often involve large accelerations of large amounts of neutron rich 
matter. Indeed the requirements for a strong source of continuous GW, at LIGO frequencies, 
places extraordinary demands on neutron rich matter. Generating GW sounds easy. Place a 
mass on a stick and shake vigorously. However to have a detectable source, one may need 
not only a large mass, but also a very strong stick. The stick is needed to help produce large 
accelerations. 

An asymmetric mass on a rapidly rotating neutron star produces a time dependent mass 
quadrupole moment that radiates gravitational waves. However, one needs a way (strong 
stick) to hold the mass up. Magnetic fields can support mountains, see for example [58]. 
However, it may require large internal magnetic fields. Furthermore, if a star also has a 
large external dipole field, electromagnetic radiation may rapidly spin the star down and 
reduce the GW radiation. 

Alternatively, mountains can be supported by the solid neutron star crust. Recently we 
performed large scale MD simulations of the strength of neutron star crust ri59l l60i . A 
strong crust can support large deformations or "mountains" on neutron stars, see also |[6T| . 
that will radiate strong GW. How large can a neutron star mountain be before it collapses 
under the extreme gravity? This depends on the strength of the crust. We performed large 
scale MD simulations of crust breaking, where a sample was strained by moving top and 
bottom layers of frozen ions in opposite directions 1,59]. These simulations involve up to 
12 million ions and explore the effects of defects, impurities, and grain boundaries on the 
breaking stress. For example, in Fig. |7]we show a polycrystalline sample involving 12 
million ions. In the upper right panel the initial system is shown, with the different colors 
indicating the eight original microcrystals that make up the sample. The other panels are 
labeled with the strain, i.e. fractional deformation, of the system. The red color indicates 
distortion of the body centered cubic crystal lattice. The system starts to break along grain 
boundaries. However the large pressure holds the microcrystals together and the system 
does not fail until large regions are deformed. 

We find that neutron star crust is very strong because the high pressure prevents the 
formation of voids or fractures and because the long range coulomb interactions insure 
many redundant "bonds" between planes of ions. Neutron star crust is the strongest material 
known, according to our simulations. The breaking stress is 10 billion times larger than that 
for steel. This is very promising for GW searches because it shows that large mountains are 
possible, and these could produce detectable signals. 

4.4. Shear viscosity of Nuclear Pasta and r-modes 

Continuous GW can also be produced by r-mode oscillations of a rotating neutron star 
1621. Consider a surface wave on a rapidly rotating neutron star that is moving slowly in a 
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Figure 7. Breaking strain (strength) of neutron star crust from MD simulations [59|. The 
shear deformation of a 12 million ion poly cry staline sample is shown. The panels are 
labeled by the strain, which is the fractional deformation. The colors in the upper left panel, 
at zero strain, show the original eight microcrystals of different orientations. The red color 
in the other three panels indicates distortion of the body centered cubic lattice. 

direction opposite to the stars rotation. This wave, in the laboratory frame, will appear to 
be moving in the direction of the star's rotation. The back reaction force on the wave from 
GW radiation will always act to slow the wave in the laboratory frame. However, in this 
case slowing in the lab frame will actually speed up the wave in the rotating star's frame and 
increase its amplitude. Thus the wave can be unstable with respect to gravitational wave 
radiation. The r-modes are collective oscillations on rotating neutron stars that can also be 
unstable to GW radiation [62] . If an r-mode is unstable the amplitude will grow large and 
rotational kinetic energy can be radiated away as GW. This will slow the rotation rate. Note, 
that the physics of large amplitude r-mode oscillations can be complicated, see for example 

162 na. 

The stability of r-modes depends on the amount of dissipation from, for example, the 
bulk and shear viscosities of neutron rich matter. If dissipation is large then the amplitude 
of the r-modes will stay small and rapid rotation of a neutron star is possible. Alternatively, 
if dissipation is small then the r-modes may be unstable and GW radiation from the modes 
may limit the rotation rate of a star. Unfortunately, the stability of r-modes has proved 
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to be a complex subject that may be sensitive to subtle dissipation properties of neutron 
rich matter, be it in a nucleon phase |l65l |66l or in more exotic quark and gluon phases 

iHIMllll. 




Figure 8. Surfaces of proton density for a pasta configuration of neutron rich matter at a 
baryon density of 0.05 fm~^. This is from a semiclassical molecular dynamics simulation 
with 100,000 nucleons 1701 . 

We give one example of a possible source of dissipation for the r-modes. The shear 
viscosity of conventional complex fluids, with large non-spherical molecules, can be orders 
of magnitude larger than that for normal fluids. This suggests that the shear viscosity of 
nuclear pasta, with complex non-spherical shapes such as the long bent rods shown in Fig. 
[8] could be large. Complex shapes arise because of coulomb frustration 171]. Here one can 
not fully satisfy both short range nuclear attraction, that correlates nucleons, and long range 
Coulomb repulsion, that anti-correlates nucleons. Pasta is expected at the base of the crust 
in a neutron star and can involve a variety of complex shapes such as flat plates ("lasagna") 
in addition to long rods ("spaghetti"). 

In ref. 11721 . we have calculated the shear viscosity of nuclear pasta using large scale 
molecular dynamics simulations. The shear viscosity is dominated by momentum carried 
by electrons and although the electron mean free path is determined by electron-pasta scat- 
tering, we find no dramatic differences from a conventional phase with spherical nuclei. 
Therefore, we find that the shear viscosity of nuclear pasta is not very different from that 
for more conventional matter with nearly spherical nuclei. However there could be other 
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sources of dissipation. We do not yet have a complete understanding of when the r-modes 
may be stable and when they are unstable. 

4.5. Dynamical Response of Nuclear Pasta 

The low energy excitation modes of nuclear pasta are important for the heat capacity, and a 
variety of transport properties, of matter deep in the inner crust. There have been some RPA 
calculations of excitation modes, see for example Ii731 1741 . Unfortunately many of these 
works use spherical Wigner Seitz cell boundary conditions that could dramatically impact 
low energy excitations. Semiclassical MD simulations allow one to calculate the dynamical 
response function, see below, and the excitation spectrum without any assumptions about 
unit cell geometries, fJ5\- 

As an example, we consider the neutrino opacity in core collapse supemovae. Neutrinos 
may scatter coherently from the pasta shapes because the shapes have sizes comparable to 
the neutrino wavelength ir76ll . The differential cross section for neutrino scattering may be 
written as follows [75l| : 

\cl{3 - x)SA{q, w) + cl{l + x)Sv{q, w)j , (6) 



dndE 47r2 

where Gp is the Fermi constant, is the neutrino energy, x = cos 9 (with 9 the scattering 
angle), and the weak vector charge of a nucleon is Ct, = —1/2 for a neutron and 
for a proton. Further, the dynamic response functions 5y(g,a;), SA{q,(^) are probed at a 
momentum transfer q and at an energy transfer uj. The axial term involving Ca = ±1.26/2 
and the dynamical spin response SA{q, w) will be discussed in a later work. Here we focus 
exclusively on the vector (density) response S{q,uj) = Sv(q,(^) that should be greatly 
enhanced by coherent effects. 

The dynamical response of the system to a density perturbation is given by 



max 



1 [T'r. 

S{q,uj) = — / S{q,t) cos{ujt) dt . (7) 



vr JO 

Here S{q, t) represents the ensemble average of the density-density correlation function 
that is computed as the following time average: 

1 1 /■T'avc 

Siq,t) = —^r^ p(q,t + s)/j(-q, s)(is . (8) 

J avc JO 

In the above expression N is the number of neutrons in the system and appropriate choices 
for Tmax and Tavo depend on the MD simulation [75 1. Finally, the one-body neutron density 
is given by p(q, t) = J2iLi exp[iq • rj(t)] with rj(t) the position of the ith neutron at time 
t. Note that because Cy^O for protons, the sum over i runs only over neutrons. Further, the 
static structure factor computed in Ref. f/Ol is easily recovered from S{q) = S{q,t = 0) = 
/o°° S{q,uj)duj. 

In Fig. [9] we show S{q,uj) for nuclear pasta at a temperature of 1 MeV 175]. The 
q = 0.05 fm^^ response shows a large peak near u = 0.8 MeV that is due to plasma 
oscillations of the charged pieces of pasta. The plasma frequency ujp depends on the charge 
and mass of the pasta, see for example Eq. [5] however Up appears to be relatively insensitive 
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w (MeV) 



Figure 9. (Color online) The dynamical response function S{q, u) of nuclear pasta versus 
excitation energy 00 =w at a density of /3 = 0.05 fm^^ and momentum transfers of g = 0.05, 
0.11, and 0.19 fm'^ Figure from ref. f!5\ . 



to the pasta shape. The peak in Fig. |9] shifts to higher excitation energies and becomes 
broader as q increases. The frequency goes up because electron screening is less effective 
at higher q, while the larger width may describe an increased coupling between Coulomb 
and nuclear excitations. 

Nuclear pasta arrises because of competition between long range Coulomb repulsion 
and short ranged nuclear attraction. Indeed at these densities there is some overlap between 
typical excitation energy scales of the Coulomb lattice and typical nuclear excitation ener- 
gies. This is reflected in the plasma oscillation peak. This peak could be important for the 
heat capacity and for transport properties. 

5. Summary, Open Questions, and Challenges for the Future 

Neutron star crust consists of neutron rich ions, electrons, and possibly a neutron gas. The 
electronic structure is simple, a degenerate relativistic gas, because of the very high Fermi 
energy. As a result, molecular dynamics simulations provide an essentially exact descrip- 
tion and can predict many crust properties. We find that diffusion is relatively fast in the 
crust because of the soft 1/r Coulomb cores between ions. This strongly suggests that 
amorphous structures will have plenty of time to crystallize, and many crystal defects will 
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diffuse away. Therefore we expect neutron star crust to be made of nearly perfect crystals 
with high electrical and thermal conductivities. Furthermore, we find that these crystals are 
very strong because of the many long range Coulomb interactions. 

There are important open questions concerning the composition of the crust. For an iso- 
lated star the composition depends on nuclear interactions including the symmetry energy. 
While for accreting stars the composition also depends on the rates of a variety of nuclear 
reactions including pycnonuclear (density driven) fusion f/Tl . 

A strong crust can support large mountains. These, on a rapidly rotating star, can gen- 
erate detectable gravitational waves. An important challenge is to understand mountain 
building mechanisms and determine what sized mountains are, in fact, present. This is 
important for many ongoing and near future searches for continuous gravitational waves. 

There are many open questions and challenges associated with nuclear pasta. Perhaps 
the most basic is "how to smell the pasta?". Which neutron star observables are sensitive to 
pasta, and which ones can demonstrate the presence of complex non-spherical pasta shapes. 
What is the shear modulus and breaking strain of pasta? What dissipation mechanisms are 
associated with oscillations of pasta? 

Finally, how do neutron stars spin so fast and what limits there rotational periods? We 
do not understand the stability of r-modes. The amplitude of r-mode oscillations is large in 
the crustal regions. Perhaps there are new sources of dissipation in the crust, such as bulk 
viscosity from an electron capture layer, that stabilize r-modes and allow rapid rotation. 
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